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ABSTRACT 

We consider a two-parameter family of cylindrical force-free equilibria, modeled to match numerical simu- 
lations of relativistic force-free jets. We study the linear stability of these equilibria, assuming a rigid impene- 
trable wall at the outer cylindrical radius /?,-. We find that equilibria in which the Lorentz factor j(R) increases 
monotonically with increasing radius R are stable. On the other hand, equilibria in which j(R) reaches a maxi- 
mum value at an intermediate radius and then declines to a smaller value jj at Rj are unstable. The most rapidly 
growing mode is an m = 1 kink instability which has a growth rate ~ (0A/~fj){c/Rj). The e-folding length of 
the equivalent convected instability is ~ 2.5 jjRj. For a typical jet with an opening angle 9j ~ few/jj, the mode 
amplitude grows weakly with increasing distance from the base of the jet, much slower than one might expect 
from a naive application of the Kruskal-Shafranov stability criterion. 

Subject headings: instabilities - MHD - galaxies: jets 



1. INTRODUCTION 

Relativistic jets in astrophysical sources have been known 
for many decades. Although their enormous power, large 
Lorentz factor and strong collimation have been widely stud- 
ied, these phenomena still lack an accepted explanation. An 
even greater mystery is the remarkable coherence and appar- 
ent stability of jets over very large length scales. This is the 
topic of the present paper. 

The most promising models of relativistic jets involve ac- 
celeration and collimation by magnetic fields with footpoints 
attached to a spinning black hole or neutron star or accretion 
disk. The forced rotation of the field lines induces a strong 
toroidal component of the ma gnetic field, which is responsible 
for accelerating the jet (e.g., iNaravan. M cKinnev. & Farmer 
20071 iTchekhovskoy. McKinney. & Nar ayan 20081 hereafter 



TMN08; and references therein). In this picture the toroidal 



component of the field dominates over other field compo- 
nents. 

Accordin g to the we ll-known Kruskal-Shafranov crite- 
rion (e.g., Bateman 1978), cylindrical magnetohydrodynamic 
(MHD) configurations in which the toroidal field dominates 
are violently unstable to the m = 1 kink instability (also called 
the screw instability). The KS criterion for instability is 
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where and B p axe, the toroidal and poloidal magnetic field 
strengths, Rj is the cylindrical radius of the jet and z is the 
length of the jet (i.e., distance from the base of the jet). Typi- 
cal jet models, including the ones described in this paper (see 
§2.2), have B^ ~ JjB p , where jj ^> 1 is the Lorentz factor of 
the jet, and they have jet angles Qj~ Rj/z ~ fsw/jj. Substi- 
tuting these scalings in equation Q), it is obvious that the KS 
instability criterion is easily satisfied in relativistic jets. We 
therefore expect astrophy sical jets to be v iolentl y unsta ble, as 
argued for example bv lBegelrrianl (1 19981) and Q (2000). Yet, 
jets in nature are apparently quite stable. How is this possible? 
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Many authors have investigated this question. They have 
used jet models with a variety of velocity profile s, geometri- 
cal sh a pes, composition and boundary condit i ons dKadomtseyl 
1966; iBatemanl [19781 iFerrari et all fl97l iBenfordl 1 198 if 
iPayne & Cohnl Il985|), and applied both analytical and nu- 
merical methods ( Istomin & Parie vl [l~996t iBegelmanl J 1998 ; 
| Lvubarskiilll999tlLill2000tlLerv et al. |2000j Appl etai|200C : 
Tomimatsu et alJl200U iMizuno et aLll2007t Moll et al.ll2008 : 
McKinney & Blandford 2008). As a result of this large body 
of work, several kinds of unstable modes have been identified: 
reflection modes, Kelvin-Helmholtz modes, current-carrying 
modes, etc. Unfortunately, it is difficult to synthesize the re- 
sults and extract universal principles. 

A fruitful approach in this field is to reduce relativistic jet 
models to their barest minimum. One such approach is to 
consider force-free jet models in which one ignores the inertia 
and pressure of the plasma and considers only charges, cur- 
rents and fields. The force-free approximation is valid when- 
ever the energy density in fields domina tes over matter en- 
ergy density, as in pulsar magnetosph eres dGoldreich & Julianl 
l!969t rRuderma n~& Sutherland!! 19751) . The force-free approx- 
imation is valid also in relativistic MHD jets, at least inside 
the fast surface (e.g., Tchekh ovskoy et aI1l2009l) . 

Theoretical studies of force-free jets have led to the iden- 
tificati on of two distinct stabi lity criteria. In a detailed anal- 
ysis, llstomin & Parievl (1 19961) snowed that cylind rical force- 
free je ts in which B z is independent of R are stable. iLyubarsklil 
( 1999) considered models with non-constant B z and showed 
that force-free jets are unstable if 



(2) 



i.e., if the poloidal field decreases with increasing dis- 
tance from the axis. We refer to equation (f2]i as 
the IPL criterion for instability. On the other hand, 
Tomimatsu, Matsuoka, & Takahashi (2001) showed via an 
approximate analysis 4 that force-free jets are unstable if 
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4 They effectively restricted their analysis to the region of the jet inside the 
light cylinder. Therefore, the flow speeds they considered were only quasi- 
relativistic at best. 
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where is the angular velocity of the field line. This criterion 
— the TMT criterion — differs from the KS criterion (HJ in 
that it explicitly accounts for rotation. It is also apparently 
very different from the IPL criterion. 

We describe in this paper a class of force-free cylindrical 
jet equilibria which clo sely matc h the numerical force-free jet 
simulations reported in TMN08. Within the context of force- 
free jets from rigidly-rotating stars, we believe that this two- 
parameter family of equilibria is generic and fairly complete. 
We study the stability properties of these equilibria and at- 
tempt to relate our results to the KS, IPL and TMT criteria 
(eqs.CD [2] and©. 

In §2 we summarize the numerical simulation results of 
ITMN08I (§2.1) and we describe an analytical force-free jet 
model which matches the simulation data very closely (§2.2). 
In §3 we carry out a linear stability analysis of these equilibria 
and show that the linear modes of the system are obtained by 
solving an eigenvalue problem involving two coupled differ- 
ential equations, with appropriate boundary conditions. In §4 
we numerically solve the equations and identify the unstable 
modes in the system. We also derive an approximate estimate 
for the growth rate of the instability. We conclude in §5 with 
a summary and discussion. We use (r, 9, (f>) for spherical coor- 
dinates and (R, <j>,z) for cylindrical coordinates. 

2. FORCE-FREE JET EQUILIBRIUM 

2.1. Structure of Force-Free Jets 

ITMN 08 considered a rigidly -rotating star of unit radius (r = 
1) surrounded by a differentially-rotating infinitely thin disk 
extending from R = 1 outward. The star was threaded by a 
uniform radial magnetic field B r , and the disk was threaded 
by a power-law distribution of vertical field, 

B z otR u - 2 . (4) 

Using a relativistic force-free co de (|Gammie et al.1 
2003V iMcKinnev & Gammid l2004t iMcKinnevI 120061 
Mignone & McKinneyl l2007t iTchekhovs koy et all 120071) . 
TMN08I numerically evolved the system and obtained the 
equilibrium configura tion of the magnetic field. 

Following TMN08, we will call the field lines that emerge 
from the star as the "jet" and the field lines from the disk as 
the "wind." The critical field line that emerges from the star- 
disk boundary defines the boundary between the jet and the 
wind. This boundary starts off at 9 = n/2 at the surface of 
the star (r = R= 1) but decreases to smaller values of 9 with 
increasing r (or z). 

We are primarily interested in the "jet" — the bundle of 
field lines attached to the central star. Since all of these field 
lines rotate at the angular velocity ft of the star, the Alfven 
surface for these lines takes the form of a cylinder — the 
"light cylinder" — with radius Ra = c/VL. Field lines become 
strongly toroidal once they are outside the Alfven surface, 
which is wher e most of the collimation and acceleration oc- 
c urs (1TMN08I) . 

TMN08 showed that the structure of the jet is strongly af- 
fected by the radial pressure profile of the region surrounding 
the jet. Specifically, if we write the radial variation of the con- 
fining pressure as p oc r~ a , the jet properties are determined 
by the value of a. In the numerical experiments, the pressure 
was caused by a force-free disk wind, and a was determined 
by the index v defined in equation (HJi according to 

a = 2(2- v). (5) 



At distance z along the axis, the cylindrical radius Rj of the 
jet is approximately given by 

Rj~ z a/4 . (6) 

For all a < 4, the jet co llimates as it moves away from the 
star (iLvnden-Belll I2006T) . As a result, at a sufficiently large 
distance from the central star (r> 1), the jet is nearly cylin- 
drical in shape. 

In the asymptotic nearly-cylindrical region of the jet, force 
bala nce in the R direction is described by the following equa- 
tion (TMNOH): 

d [B 2 -E 2 \ (Bl-E 2 \ (b 2 -E 2 \ 

where B is the total magnetic field strength, B^ is the toroidal 
field strength, and B p is the poloidal field strength: 

B =^ B l + B l B p = yj# R +& v (8) 
The electric field E is given by 

E = -(SlR/c)4>xB, (9) 

where SI is the angular velocity of the field line. The electric 
field has only a poloidal component: E p = VLRB p /c. 

Each of the three terms on the left-hand side of equa- 
tion (O represents a force in the -R direction. The quantity 
(B 2 -E 2 )/8ir is the pressure of the force-free fluid in the co- 
moving frame; therefore, the first term describes the inward 
force due to the gradient of pressure. The second term arises 
from the toroidal curvature of the field line. The toroidal 
magnetic field B^ contributes an inward force due to "hoop 
stress," while the poloidal electric field E contributes an out- 
ward force. 5 The third term in (jTj gives analogous contribu- 
tions from the poloidal curvature of the field line, where R c 
is the poloidal radius of curvature; once again there is an in- 
ward force due to the poloidal magnetic hoop stress and an 
outward force due to the electric field. Note that the contribu- 
tions involving E are important only for relativistic flows. In 
standard non-relativistic MHD, one neglects these terms and 
k eeps only the terms involving B. 

TMN08 showed that models with a > 2, i.e., v < 1, are 
good analogs of relativistic jets found in nature (especially 
the jets of gamma-ray burst). Figure |2~T1 shows numerical re- 
sults cor responding to the "fiducial" force-free simulation in 
TMNOl with v = 0.75 (equivalent to a = 2.5). Panels (a) and 
(b) show results for a relatively near region of the jet at z = 10 2 , 
and panels (c) and (d) show results for a more distant region 
at z = 10 7 . In each case, the abscissa corresponds to the cylin- 
drical radius R normalized by the local "jet radius" Rj, which 
is the cylindrical radius of the last jet field line that separates 
the jet from the surrounding disk wind. 

The main features of the numerical solution are as follows. 
First, from panels (a) and (c) we see that B p is essentially con- 
stant inside the jet, showing hardly any variation with R. As 
we show in AppendixlAl this is required for force-free jet so- 
lutions that smoothly connect to the central compact object. 6 

5 By equation |9}, E is directly proportional to f2, so the latter term results 
from rotation and may loosely be viewed as a sort of "centrifugal force" (V. 
Beskin, private communication). 

6 Asymptotic force- free jet configurations w ith non-constant profiles of B p 
are certainly possible l lstomin & Pariev 1996), but there exists no solution 
that would smoothly connect them to the compact object. 
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FIG. 1 . — Numerical results from a force-free jet simulation with a = 2.5 ("fiducial model" corresponding to v = 0.75 inlTMN08). [Panel (a)]: Field components 
as functions of normalized cylindrical radius R/Rj at z = 10 2 . The solid horizontal line shows B 2 (R)/B 2 (0), the short-dashed line shows [B 2 i (R)—E 2 (R)]/B p (0) 
and the long-dashed line shows E 2 (R)/B 2 (0). The dotted lines are the corresponding results for the analytical model with 7,,, = 6, R m = 1.2 (Model A, §2.2.2). 
The vertical solid line shows the boundary between the jet and the external confining medium. [Panel (b)]: Lorentz factor 7(R) (solid line), and the two 
approximations, 7i(K) (short-dashed line) and "ji(R) (long-dashed line), at the same z. The dotted lines are from the analytical model. [Panels (c), (d)]: Similar 
to (a), (b), but at z = 10 7 . The dotted lines in these panels correspond to the analytical model with 7,„ = 2600, R m = 0.18 (Model C, §2.2.2). 



Second, since VL is constant and E oc VlRB p (eq. we have 
E 2 ocR 2 . Third, is almost equal to E and so also varies 
primarily as R 2 . Fourth, \B^\ is slightly larger than E with 
B^ — E 2 oc R 4 . Because of this property, the first two terms 
in equation O both give an inward force. Therefore, we ob- 
tain the fifth feature of the solution, viz., the third term in 0, 
which involves the poloidal curvature of the field line, is im- 
portant for force balance. This is the only outward force in 
the balance equation - it is outward because E is of order B^ 
and is much greater than B p outside the light cylinder. This 
force has to balance the other two terms. 

The velocity of a force-free flow is usually identified with 
the drift velocity, 

V ExB 

-c = ir> (10) 

and the Lorentz factor is defined correspondingly. Since E ■ 
B = 0, it is easily shown that 



B 1 



7 



B 2 -E 2 ' 



(11) 



Panels (b) and (d) show the variation of 7 as a function of the 
normalized cylindrical radius R/Rj in the numerical model. 



TMN08 derived two approximate relations for 7, 

7 i = [l + W C ) 2 ] 1 / 2 , 
l2 = OR c /R) 1 ^ 2 , 



(12) 
(13) 



and they showed that the net 7 of the fluid is given to good 
accuracy by the following simple formula, 



1 _ 1 1 

T if 72 



(14) 



Along each field line, 7 is initially determined mainly by ro- 
tation, and_so_7 « 71 . This is a region of efficient acceleration 
which TMN08 called theirs? acceleration regime. However, 
beyond a certain distance from the star, the effect of poloidal 
field line curvature becomes important, and 7 switches to the 
less efficient 72, the second acceleration regime. 

Relatively near the star, all field lines in the jet are in the first 
acceleration regime and j(R) behaves like 71 (eq. fl2l) . as seen 
in panel (b). Specifically, 7 increases more or less linearly 
with R and reaches its maximum value at the edge of the jet at 
R = Rj. However, when we consider the jet at a larger distance 
from the star, some of the field lines have already switched 
to the second acceleration regime, where 7 ~ 72 oc l/R (see 
eq. [13] coupled with eq. l2"Tlbelow). We then have the results 
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shown in panel (d), where the maximum Lorentz factor occurs 
at some radius R„, inside the jet, not at the boundary; we have 
7 ~ 7i o c RforR< R m and 7 ~ 72 oc R~ l for R„, <R< Rj. 
TMN08 discuss in detail the physics of the two acceleration 



regimes. 



2.2. Analytical Jet Model 



The axisymmetric numerical jets models described in §2.1 
have magnetic and electric field components that are functions 
of both R and z. This is not convenient for linear perturbation 
analysis. Since the numerical models are nearly cylindrical at 
large distance (i.e., dR/dz^. 1), we consider now an ideal- 
ized jet equilibrium model which is perfectly cylindrical and 
in which all quantities are functions only of R. We choose the 
following specific functional forms: 



Bqr = 0, 



B ^ = - [2( 7 „ 2 ,- \){R/R m f + (R/R m ) A ] 1/2 : 



B fe = exp [-3 J R 2 /4( 7 „ 2 , 

E 0R =- [2(7^-1)] 1/2 
£b</> = 0, 
Eo, = 0, 

R c = 2( n l-\)R 2 J2>R. 



l)R 2 m ] =g(R), 
(R/R m ) = -h(R), 



(15) 

-f(R), (16) 
(17) 

(18) 
(19) 



(20) 
(21) 

The zeros in the subscripts are meant to indicate that all these 
quantities refer to the unperturbed model. The model has two 
parameters, -f m and R m , whose meanings are explained below. 
For simplicity, we have chosen units such that Bq z = 1 at the 
jet axis (R = 0). Also, we have assumed that Bo z and Q are 
positive, so both Bq^ and Eqr are negative, i.e., magnetic field 
lines are swept backward with respect to the rotation and the 
electric field is pointed radially inward. With this choice of 
signs, the three functions, f(R), g(R) and h(R), are positive. 
Note that, in all cases of interest, g(R) is practically equal to 
unity. The particular exponential form given in equation ( fTvT > 
is designed to handle small higher-order terms in the force 
balance equation (1221 , but the deviations of g(R) from unity 
are tiny and unimportant. 7 

By direct substitution it is easily verified that the above 
model satisfies the radial force balance equation Q. Under 
cylindrical symmetry, this equation takes the form: 

d_ ( B U + B lz ~ E ^ 
dR \ 



8tt 
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4ttR 
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(22) 

The first two terms are positive, i.e., both represent inward 
forces, with the first term providing twice as much force as 
the second. The third term is negative and its magnitude is 
equal to the sum of the other two terms. 

An important feature of the above model is that the out- 
ward force from the third term involves the poloidal curva- 
ture radius R c . Technically, a perfectly cylindrical model has 
R c — > 00. To get around this problem, we treat R c as an ex- 
ternally imposed property of the solution which is adjusted 
so as to reproduce the poloidal curvature force present in the 
numerical jet model. In other words, even though we have 
straightened out field lines in the z-direction by enforcing 
cylindrical geometry, we still retain the effect of poloidal cur- 
vature by means of an artificial external force. This procedure 

7 As a test, in the stability analysis described later we h ave done the cal- 
culations both with the full expression for g(R) given in eq. i\7\ and with the 
simpler choice g(R) = 1 . The results are practically the same. 



is analogous to the widely-us ed shearing sheet approxima- 
tion in accretion di sk studies (Goldre ich & Tremaindll978l : 
iNaravan et al.ll 1987b in which fluid streamlines are straight- 
ened out in the azimuthal direction, but the effect of az- 
imuthal curvature is still retained via a Coriolis force. Note 
that, apart from the extra term due to poloidal field curva- 
ture, equation (j7]) is identical to the standard balance condi- 
tion d erived in other papers i n the literature, e. g., equation 
(6) in Istomin & Parievl (1 19961) or equation (16) in lLyubarskiil 
(1999). 

To get a better idea of the nature of the above analytical 
model, we now make a couple of simplifications. As already 
mentioned, Bq z is practically independent of R inside the jet. 
Also, for highly relativistic jets, we invariably have B^-E 2 <C 

B?,. We therefore replace equations ( TTSI l. ( flTl and ( fT8l by the 
following simpler formulae 



Bo^kEqh 



[2(7 2 -l)] 1/2 (fl/fl ffl ) 
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(R/R,n) 4 



(23) 
(24) 
(25) 



By equation (O, the angular velocity of rotation of the field 
lines is given by 



= -cE OR /B 0z R « J2( 7 2 -\){c/R m ), 



(26) 



and is the same for all field lines, as required for a rigidly ro- 
tating star at the base of the jet. 8 Using these simpler expres- 
sions, we obtain the following result for the Lorentz factor: 

1 B 2 -E 2 B 2 +B 2 Z -E 2 R 



1 2 (R) 



B l</, + B 0z 



\ + (R/R m f 



\+2{ 1 2 m -l)(R/R m ) 2 



l + (flR/c) 2 3R C ' 



(27) 
(28) 



where we have made use of equations (l2Tb and d26i >. We thus 
reproduce the result given earlier in equation (|14b . 

It is easily shown that 7 reaches a maximum at R = R„, and 
that its value at this radius is equal to 7,,,. Thus, the two model 
parameters R m and 7,,, allow us to control the basic features of 
the equilibrium. 

2.2. 1 . R,„ > Rj: Maximum Lorentz Factor Located at the Jet 
Boundary 

A jet in which all field lines are in the first acceleration 
regime has its maximum Lorentz factor at the boundary of 
the jet, R = Rj. This corresponds to choosing R m > Rj in the 
analytical model, so that the term (R/R m ) A in the numerator of 
equation d27T ) can be neglected. In this case, the profile of 7 
has two segments: the region of the jet inside the light cylin- 
der (R < Ra = c/Vl) which is not accelerated very much, and 
the region outside the light cylinder which has 7 increasing 
linearly with radius, 



7 (i?)( 



I 1 ' 
{R/Ra 



(R/RjYfj, 



R<R A = R„,/y/2(ri-l), 
R A <R< Rj, 



8 Since we have designed our analytical model to match the numerical 
models of TMN08, all of our models have constant Q(R). It would be 
straightforward to generalize the model to non-constant Q(R), using addi- 
tional parameters. 
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FIG. 2. — Profiles of 7 vs R/Rj for the analytic Models A, B and C. 

where is the Lorentz factor at the jet boundary, 

jj « 7l (/? ; ) = [l +2( 7 2 - 1)0R ; -/X) 2 ] 1/2 . (29) 

The dotted lines in panels (a) and (b) in Fig. 12. II show the 
dependences of various quantities as a function of R/ Rj for a 
model with 7,, = 6 and R m = l.2Rj. The agreement with the 
numerical simulation results at z = 10 2 is striking. We call the 
analytical model with this particular choice of 7,, and R m as 
Model A: 

Model A : 7m = 6, R m = 1 2Rj. (30) 

The solid line in Fig. |2]shows the variation of 7 as a function 
of R for this model. 

2.2.2. R„, < Rj: Maximum Lorentz Factor Located Inside the Jet 

As we described in §2.1, ajet in which some field lines have 
switched to the second acceleration regime has its maximum 
Lorentz factor inside the jet. This means R„, < Rj. In this case, 
the Lorentz factor jj at the jet boundary is roughly equal to 

7j « 72(^0 = [2(7* - D] 1/2 RjRj- (3D 
Now the profile of 7 has three segments: 

fl, R<Ra, 
j(R) « I R/R A w (R/R m )l,n, R A <R< R m , (32) 
[ (R m /R) lm « {Rj/Kftj, R m <R< Rj. 

The dotted lines in panels (c) and (d) of Fig. 12.11 show 
model results corresponding to 7,,, = 2600 and R m = OARRj. 
We find excellent agreement with the numerical results at 
Z = 10 7 . We call the analytical model with these values of 
7m and R m as Model C. For completeness we also consider 
a less extreme model called Model B in which 7m = 6 and 
Rm = 03Rf 

Model B: 7 m = 6, R m = 03Rj, (33) 
Model C: 7,, = 2600, R m = 0.18Rj. (34) 

The dashed and dotted lines in Fig. [2] show the variations of 7 
as a function of R for these two models. 



3. LINEAR PERTURBATION ANALYSIS 

We now consider linear perturbations of the cylindrical 
equilibrium described in §2.2. The unperturbed state has mag- 
netic and electric fields 



Bo = BqrR +Bo4><fr+Bo z z, 
Eo = EqrR +Eo,j > <p+Eo z z, 



(35) 
(36) 



where the various components are given by the expressions 
in equations <fT3T> — (f20b. As mentioned previously, we choose 
Bo z and f2 to be positive, so Bq^ and £bs are negative. The 
unperturbed current and electric charge are 
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Note that we have included terms involving R c in the un- 
perturbed current and charge density. These terms describe 
the contributions of poloidal field curvature to the quantities 
V x Bo and V • Eq, respectively. By including these terms, 
we retain the forces associated with poloidal field curvature 
without actually having curved field lines in the model. 

3.1. The Eigenvalue Problem 
We now consider small perturbations. Let us write the per- 
turbed electric field as E = Eq+Ei, where E\ is a small pertur- 
bation of the form 

Ei = [E lR (R)R+E l4 ,(R)4>+EUR)z] exp(-iut + im(f>+ikz). 

(37) 

Let us similarly write B = Bq+B\, J = Jq+J\, p = po + p\. Each 
of these small perturbations can be expressed in terms of the 
perturbed electric field via Maxwell's equations. From 



IdB - 

c at 



we obtain 



From 



we obtain 



Finally, from 



we obtain 



•Vx& 



1 dE - 4tt - 
-— = VxB J, 

c at c 



IUJ -> IC -> -> -> 

-Ei-- — Vx(Vx£i). 



4tt 4irui 



V-E =47rp, 



Pl = -V.£, 



(38) 



(39) 



(40) 



(41) 



(42) 



(43) 



Since the perturbed system is force-free, it must satisfy E ■ 
B = 0. The zeroth order terms satisfy this trivially (as they 
should). The first order terms give the condition E\ -Bq + Eq ■ 
B\ = 0. Substituting for the various quantities, this condition 
allows us to solve for E\a, in terms of E\ z : 



E\d> = C\E\ 



(44) 
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where the function C\ is given by 

tuRg — cmh 
Cl ~ uRf-ckRh' 
Successive differentiations give 

Ei < t > = C\E Xz + C x E u , 

II II II II 

^i^ = ^'i^'iz + 2C 1 £' lz + C 1 E\ z 



(45) 

(46) 
(47) 



We now consider the force balance condition: pE + ( 1 /c)J x 
B = 0. The zeroth order terms give 



(48) 



PoE +-JqxBo = Q, 
c 

which is simply the equilibrium force balance condition i 
Notice that the poloidal curvature terms in Jo and po are neces- 
sary to satisfy equilibrium in the unperturbed solution. From 
the first order terms in the force balance equation we obtain 

PiEo + poEi + -fx x B +-J x B x = 0. (49) 

c c 

The (f> component of this equation gives a relation between 
E^, E\$, E lz and E\ z . Eliminating E i( , using equation (|46| >. 
we obtain a first-order differential equation for Ei z (R): 

D 1 e' 1z + D 2 E 1z + D 3 E ir = 0, (50) 

where D\, D 2 and D3 are functions of R. The expressions 
are relatively long and we give them in Appendix [B] The z 
component of d49l has no new information; it just gives back 
equation ( f44b . The R component, however, gives a new rela- 
tion between the various components of E\ and their deriva- 
tives. Eliminating E l( p, E 1( p, E lz and E lz using equations (|46*T i. 
( |471 l, d50l > and the derivative of ( [50) , we obtain a differential 
equation for E lR (R): 



D 4 E lR + D 5 E lz + D 6 E lR = 0, 



(51) 



where D4, D5 and Df, are again functions of R and are given 
in Appendix IE1 

We have thus reduced the linear mode analysis problem to a 
pair of first-order differential equations, d50l >, ( Bil l, for Ei z (R) 
and Eir(R). For convenience, we write down the two equa- 
tions again: 

D 2 D 3 



Eu Eir, 

Di K Di ' 



~D 4 



-Eu- 



D 4 



(52) 



(53) 



These equations constitute an eigenvalue problem, where u> 
is the eigenvalue. By numerically solving the equations with 
appropriate boundary conditions, we obtain u> for given values 
of m and k. 

The singular points of the equations are located at the radii 
where D\(R) and D^iR) vanish. Anticipating later discussion, 
we write down here the expression for the quantity D1D4: 

" (uRf- ckRhf + (ojRg - mch) 2 - (ckRg - mcfj 1 ' 



0^4 = -— 

ujR 



(ujRf-ckRh) 
Also, from equation ( fTOb , the perturbed velocity is 

v\ Ei x Bo + Eo x B\ 
~ = ~B 2 ' 



(54) 



(55) 



which in component form gives 

v\r _ gEjf+fEiz 
c 

Vl<t> 

c 

vu 
c 



f 2 + 8 2 ' 
gE iR -hB l: 

f 2 +g 2 
fE lR + hB lq 



(56) 
(57) 
(58) 



f 2 + 8 2 ^ 

We now consider boundary conditions. A physically valid 
perturbation will be well-behaved on the axis (R = 0) and 
will satisfy suitable boundary conditions at the jet boundary 
(R = Rj). The condition on the axis is different for axisym- 
metric (m = 0) and non-axisymmetric (\m\ > 1) perturbations, 
so we consider each of these cases in turn. At the jet bound- 
ary, we assume that the jet is constrained by a "rigid wall" and 
we write down the corresponding boundary condition. In the 
following, we employ the specific forms of f(R), g(R), h(R) 
given in equations ( TT6l >, (fl7l l. ( TT8l . 

3.2. Boundary Condition on the Axis: m = 

Setting m = and substituting the expressions for f(R), 
g(R), h(R) in D\ -D§, we find that the leading terms of the 
differential equations ( |52l , ( 1331 ) at small R are given by 



where 



E\z~ -rrE\ z + a zR Ei R , 

K 

/ a R , a RR 



2uj 
ck ' 
i(c 2 k 2 -uj 2 ) 

Aiu 
k(ck-Lu) ' 
(ck+2uj) 



a zR -- 



a Rz -- 



a RR =- 



ck 



(59) 
(60) 

(61) 

(62) 

(63) 

(64) 
immedi- 

(65) 
(66) 



Requiring the perturbation to be analytic as R 
ately gives the following solution near the axis, 

E lz = KR 2 , 

He 

E\ R = —K R : 

(ck—uj) 

where K is an arbitrary normalization constant. 

3.3. Boundary Condition on the Axis: \m\ > 
When m ^0, we obtain 

a zz = 0, (67) 
itnA(ck— oj) 



a zR -- 



(Acm — ujR m ) ' 
im(Acm — LoR m ) 



a Rz — - 

A(ck-uj) 

a RR = -l, 

where the constant A is defined to be 

2 , vi 1/2 



(68) 

(69) 
(70) 



A= [2( 7 2-l)] (71) 
The physically relevant solution close to the axis is then 

E lz =KRW, (72) 



£ifi = -^Tsgn(m) 



i(Acm — ujR m ) 
A(ck-uj) 



m I - 1 



(73) 
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3.4. Boundary Condition at the Jet Boundary: Rigid Wall 

We assume that our cylindrical jet is terminated at R = Rj 
by a rigid impenetrable wall. By impenetrable we mean that 
no energy flows across this boundary, either out of or into the 
jet, i.e., the Poynting flux lies in the <f>— z plane. Equivalently, 
the velocity vector has no radial component. 

The equilibrium Poynting flux of course lies in the <f> — z 
plane. The perturbed Poynting flux is proportional to E\ x 
Bq + Eq x B\. Since Eq is parallel to R, the term Eq x B\ is 
automatically in the <fi— z plane. The term E\ x Bq will also 
be in this plane if E\ is precisely radial, i.e., both and E\ z 
vanish. By equation d44i >. is proportional to E\ z . We thus 
obtain the following boundary condition at the outer wall: 



R = R, 



(74) 



4. NUMERICAL RESULTS 



We have computed frequencies of modes by numerically 
solving the differential equations (l52l and (1531 . along with the 
boundary conditions described in §§3.2-3.4. For each choice 
of k and m, a countable infinity of solutions exists which may 
be ordered by the number of zeros of Re[Ei z (R)], not count- 
ing zeros at the boundaries. 9 The lowest-order solution (the 
"fundamental mode") is such that Re[£ , i ; (/?)] has no zeros be- 
tween R = and R = Rj, the next solution has one zero inside 
the jet, and so on. In the following we identify each mode by 
its radial mode number n which we define to be the number 
of zeros. As one might expect, the mode frequency increases 
with increasing n. 

We solve for the frequencies via a shooting method. We 
start with a guess value of u>, make use of the expressions 
given in §3.2 or §3.3 (depending on the value of m) to set up 
the initial solution at small R, and integrate equations 02] ) and 
d53l to R = Rj. We then adjust uj in the complex plane until the 
outer boundary condition given in §3.4 is satisfied. The only 
subtle point is that the quantities D\ and D\ appear in the de- 
nominators of various coefficients in equations d52l and ( |53l , 
and so their zeros correspond to poles in the solution. To avoid 
these poles, we treat R as a complex variable and integrate 
the equations over a "safe" trajectory in complex-/? space. 
Since the solution is analytic, the exact track that we follow 
is unimpor tant so long as i t lies a bove all singularities in the 
/?-plane. Istomin & Pariev ( 1996) give a detailed discussion 
of this topic in connection with current-driven instabilities in 
force-free jets. The reader is also referred to standard discus- 
sions of this point i n plasma physics texts in the context of 
Landau damping, or Goldreich, Goodm an. & Naravanl(ll986h 
for a discussion in the context of accretion disk instabilities. 

4.1. Axisymmetric Modes: tn = 

Figure [3] shows the dispersion relation — the variation of u> 
with k — of axisymmetric modes (m = 0). Results are shown 
for Models A, B and C (eqs. 1501 15311311 for three radial mode 
numbers: n = 0, 1, 2. For all k, we find that the mode fre- 
quency is real, which means that all these modes are stable. 

For large values of kRj, the mode frequency asymptotes 
to uj = ±ck, so the modes behave like electromagnetic waves 
moving parallel or anti-parallel to the z-axis. At small k, how- 
ever, the frequency asymptotes to a constant value. There is 
thus a minimum frequency for propagating modes inside the 
jet. The minimum frequency is of order the inverse of the 

9 Re() stands for the real part of a complex quantity. 
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FIG. 3. — Dispersion relation for axisymmetric modes (hi = 0) in Model A 
(solid lines), Model B (dashed lines) and Model C (dotted lines). From below 
the curves correspond to radial mode numbers n = 0, 1,2. All the modes are 
stable. 
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FIG. 4. — Imaginary part of ui for modes in Model B with m — 1, n = 
Growing modes have lm(u)) > 0, while decaying modes have lm(uj) < 0. 



light-crossing time across a radial wavelength of the mode 
(e.g., uj m \ n ~ 2-Kc/Rj for the mode with n = 0). 

We find that the dispersion relations of modes with positive 
and negative k are not quite the same. The difference arises 
because the background has a non-zero velocity in the z direc- 
tion, which breaks the symmetry between waves propagating 
towards +z and — z. The effect is, however, quite weak. 

4.2. Non-Axisymmetric Modes: m = ±1 

The most interesting modes are those with m= ±1. These 
modes are stable in Model A, but unstable in Models B and 
C. 
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FIG. 5. — Eigenfunctions coiTesponding to growing modes in Model B with 
m = 1, n = and kRj = 5, 10, 15, 20 and 25. The real part of E\ z is plotted. 

Figure[4]shows Im(aj) 10 as a function of kRj for a sequence 
of modes in Model B; the modes correspond to m = 1, n = 0. In 
this sequence, modes with < 0.65 and those with kRj > 28 
are stable and have uj real. However, for 0.65 < kRj < 28, we 
find a pair of modes with complex values of uj. The branch 
with Im(oj) > corresponds to growing modes, and the branch 
with Im(cj) < to decaying modes. 11 

Figure [5] shows eigenfunctions corresponding to a few of 
the growing modes. Plotted is Re(£ , i < .) as a function of the 
scaled radius R/Rj- The mode corresponding to kRj = 5 is 
representative of all modes with kRj < 5. These modes have 
eigenfunctions with no zero crossings between R = and R = 
Rj. By our definition, the modes correspond to n = 0. Each 
of the remaining eigenmodes in Fig. |5]has a pronounced dip 
in Re(Ei z ) which causes a zero crossing. These dips result 
from a singularity in the equations, as we discuss below. If 
we discount the singularity-induced zero crossings, then these 
eigenfunctions may also be identified as n = modes. 

Figures [6] and [7] show similar results for Model C. Grow- 
ing modes (and their decaying counterparts) are present for 
m = 1 and all kRj < 2.1 x 10 4 . Modes with m = — 1 are also 
unstable (see Fig. |6]l 12 . Figure [7] shows a few eigenfunctions. 
All modes with kRj < 1.3 X 10 3 have eigenfunctions with the 
standard n = shape (see the mode with kRj = 1.25 x 10 3 in 
Fig. |7). For larger values of k, the eigenfunctions develop 
negative spikes due to the presence of a singularity (see Fig. 
[7J. However, we still view them as n = modes. 

We have determined numerically that the singularities 
which cause the dips in the eigenfunctions are due to zeros 
in the function D^R) defined in §3.1. This function appears 
in the denominator of the differential equation (l53l . and hence 
its zeros behave like poles. 13 

Equation d54l i gives the analytic form of the quantity D\D\. 

10 Im() refers to the imaginary part of a complex quantity. 
1 ' We discuss the sudden jump in the value of Im(aj) for the decaying 
branch at the end of §4.2. 

12 In the case of Model B, modes with m = — 1 appear to be stable, and only 
the m = +1 modes show an instability. 

13 In c ontr ast, although the function D\(R) appears in the denominator of 
equation {52}, its zeros do not cause a real singularity since the terms Di and 
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FIG. 6. — Imaginary part of lu for growing modes in Model C. The solid 
line corresponds to modes with m = 1, n = and the dotted line corresponds 
to modes with m = — 1, n = 0, 
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FIG. 7. — Eigenfunctions corresponding to growing modes in Model C with 
m=l, n = QmdkRj = 1.25 x 10 3 , 2.5 x 10 3 , 5 x 10 3 , 10 4 and 2 x 10 4 . The 
real part of E lz is plotted. 

Since the modes of interest to us have Re(w) very nearly equal 
to ck, let us substitute u> = ck in this equation. Then, setting 
D1D4 equal to zero gives the following relation between the 
wavenumber k of the mode and the radius R sm „ of the singu- 
larity: 

kR _ w[/(flsing) + Musing)] 

SmS " g(fismg) + sgn(m) y/g 2 (R™ S ) + [/^sing) - h 2 (R smg )] ' 

(75) 

Figure [8] shows the position of the singularity R smg as a func- 
tion of kRj for modes with m = ±1 in Models B and C, as 
calculated with this equation. For comparison, the dots show 
the radii at which the functions Re [£1 ;.(/?)] reach their minima 

Z>3 also go to zero at the same locations. 
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FIG. 8. — Location of the si ngula rity K s j ng as a function of kRi in Models 
B and C, calculated using eq. )75t . Solid curves correspond to modes with 
m= 1, n= and dotted lines to modes with m = —i, n = 0. Solid dots show 
the locations of minima in the eigenfunctions plotted in Figs.|3]and|2] 

in the eigenfunctions plotted in Figs. |5]and|7] The agreement 
between the analytical curve and the dots is excellent, show- 
ing that equation (l75l > captures the physics of the singularity. 

From Fig. [8] we see that the singularity lies inside the jet 
Casing < Rj) only for a finite range of k above a certain min- 
imum value. For values of k smaller than this minimum, the 
singularity is outside the jet (for very small k it is well out- 
side the jet). In the case of Model B, the singularity enters the 
jet from outside when kRj ~ 5 and it disappears (for m = +1) 
at R = when kRj ~ 28. This is the primary range of k over 
which an unstable mode is present. At kRj ~ 28, the singular- 
ity is barely present near the center of the jet and we have a 
very weakly growing mode. With decreasing k, the singular- 
ity moves outward and the growth rate of the mode increases 
(Fig. |4j. At kRj ~ 5, when the singularity reaches the wall, 
the growth rate is close to its maximum value. At yet smaller 
values of k, the singularity moves outside the outer wall, but 
its presence is still felt and there is continued instability. The 
growth rate however decreases with decreasing k. 

A similar pattern is seen in Model C. Unstable modes are 
present only for kRj < 2 x 10 4 . With decreasing k the growth 
rate increases and reaches its maximum approximately when 
the singularity reaches the jet boundary (/? s i n£ , = Rj), which 
happens at kRj ~ 1.3 x 10 3 . In contrast to Model B, however, 
the growth rate remains large even for smaller values of k, and 
the instability survives down to k — > 0. 

We finally discuss the peculiar behavior of Im(cj) in the de- 
caying branch of modes in Fig. [4] As we mentioned earlier, 
in numerically solving for the eigenvalue we must integrate 
the differential equations (l52l and d53l along a path in the 
complex-/? plane that lies above the poles in the solution. For 
growing modes, the pole is located below the real R-axis. We 
can therefore integrate along the real R-axis without any dif- 
ficulty. For decaying modes, the pole is above the real R-axis 
and now we must choose the integration path with care. If the 
singularity has Re(/? s i„ g ) > Rj, i.e., the singularity is outside 
the jet, there is no problem and we can simply integrate along 



the real axis. However, when < Re(/? s j ng ) < Rj, we have to 
deform the integration path. In our calculations, we integrate 
from R = along a path with lm(R) = Re{R) until the point 
lm(R) = Re(R) = Rj and we then integrate down to R = Rj. 
The jump in Im(cj) in Fig. |4]is the result of the singularity 
moving into the jet. To the left of the break, the singularity is 
located at R > Rj. Here the eigenvalues of the growing and 
decaying modes are complex conjugates of each other. How- 
ever, to the right of the break, the singularity has moved inside 
the jet (R < Rj) and now the complex conjugate symmetry is 
broken. 

We note that eigenfunctions and eigenvalues of decay- 
ing modes are not very meaningful. This can be shown 
from an initial condition analysis along the lines of Lan- 
da u's treatment of plasma damping. The reader is referred 
to llstomin & Parievl (1 19961) for a detailed discussion of this 
topic. 

4.3. Why m = ±1 is Special 

We have not exhaustively explored modes with \m\ > 1. 
However, in spot tests with various choices of m, n and kRj in 
Models A, B and C, all modes were found to be stable. We 
believe that, if at all, there are only weakly unstable modes 
for |m| > 1; there is no sign of the kind of vigorous insta- 
bility described in §4.2 for modes with m = ±1. So why is 
m= ±1 special? The answer to this questio n is well-know n 
in the magnetic confinement literature (e.g., Bateman 1978b . 
We discuss it briefly here for completeness. 

Consider the radial component of the perturbed velocity v^? 
near the axis of the jet. Equation (TS6b gives the expression for 
v\r in terms of the perturbed electric field components E\ z and 
Ei,p, and equation d44l) shows the relation between these two 
field components. For small values of R near the axis, we have 



g(R)nl, f(R)^h(R)=A^-, 
where the quantity A is defined in equation d7ll . and 



C, 



(ujR m — Acm) 
AR(uj-ck) 



(76) 



(77) 



Consider first modes with m = 0. Equation d65l ) shows that 
E lz w KR 2 near the axis. Substituting this in equation d56l ) and 
using the other approximations given above, we find 



ujR„ 



A(oj-ck) 



KR + 0(R 3 ). 



(78) 



By symmetry, the velocity goes to zero on the axis, and the 
flow consists of a simple radial divergence. 

Consider next modes with m ^ 0. Using equation d72b for 
E\ z , we find 



v« « (UJR "' Acm) KRW- l cosmq> + 0(RW +l ), 
A(lu - ck) 



(79) 



where we have included cos m<p to show the angular depen- 
dence of the mode. The leading term goes like /?M -1 , which 
corresponds to R° when m = ±1. This means that the mode 
has a finite radial velocity, and hence a finite radial displace- 
ment, on the axis when \m\ = 1. The cos0 dependence of vir, 
coupled with the fact that v 1( p has the same amplitude but a 
sine/) dependence, ensures that the velocity vector is unique 
and analytic at R = 0. By writing the velocity vector in carte- 
sian coordinates, it is easily seen that the complex phase of 
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vm determines the orientation of the velocity vector in the xy- 
plane. If we consider values of |m| > 2, the velocity vanishes 
on the axis, just as in the case of m = 0. 

This then reveals what is special about \m\ = 1 modes. 
These are the only modes in which fluid perturbations com- 
municate across the axis and cause the jet to shift bodily 
across the axis. In modes with \m\ = 1 the center of mass of the 
jet itself shifts into a spiral shape, which is the characteristic 
feature of the kink or screw mode. For all other values of \m\, 
the center of mass remains on the axis and the perturbations 
are concentrated on the outside. 

In helical MHD configurations in the laboratory, the \m\ = 
1 kink mode is known to be highly unstab le and to be th e 
greatest threat to the stability of equilibria (Bat emar3ll978l) . 
Not surprisingly we see the same feature in our force-free jet 
equilibria. 

4.4. Growth Rate of the Instability 

The growth rate of the fastest growing mode is a matter of 
practical interest since it limits the lifetime of an unstable sys- 
tem. As discussed in §4.2, for the models we have considered 
here, the most unstable mode generally has a singularity close 
to the outer wall: R s \ ng ~ Rj. Knowing this, we estimate here 
the fastest growth rate by assuming that the pole is located at 
^sing = l.lRj. (We locate the singularity slightly outside the 
jet, since this speeds up the numerical integrations consider- 
ably.) 

Given an assumed value of R s \ ng , we can substitute this 
value in equation ( TTSb and make use of the expressions for 
f(R), g(R), h(R) given in §2.2. Recalling that Models B 
and C are in the regime described in §2.2.2, we note that 
f 2 (Rj)-h 2 (Rj) 3> g 2 (Rj)- In addition, / and h are nearly equal 
to each other and jj is given by equation (|3T1 >. We then find 
that kRj pa 1 .6jj. Also, the real part of the frequency is nearly 
equal to ck. Thus, we estimate 

I.67; 1.67,c 
Mode with R sin „ = 1 . IRj : k pa , Re(u) pa 'J- . 

R j R j 

(80) 

These estimates should apply to the fastest-growing mode. 
The mode with the maximum growth rate in Model B has 
Re(w) pa kRj = 4.4. Since Model B has = 2.5, equation 
( T80T > predicts kRj pa 4.0, which is close. Similarly, the mode 
with the maximum growth rate in Model C has Re(w) pa 
kRj = 1000, whereas equation (l8TT l with jj = 660 predicts 
kRj pa 1060. We see that the approximate formula (l80l is quite 
good. 

Our numerical results indicate that the growth rate Im(cj) 
of the fastest growing mode is proportional to Re(o;)/7?. 
We also know that unstable modes are present only when 
R,„ < Rj\ for instance, Model A with R m = 1 .2Rj has no unsta- 
ble modes, whereas Model B with R m = 0.3/?; and Model C 
with R„, = 0.18/?j both have unstable modes. With these clues 
in mind, we obtain the following empirical estimate for the 
growth rate of the fastest-growing mode: 

0.4 / 2R m \ c 
Mode with R sin „ = 1 . 1R, : Im(w) pa — 1 — . 

n \ Rj J Rj 

(81) 

The coefficients 0.4 and 2 are very approximate (to emphasize 
this, we give only the leading digit for each). Nevertheless, as 
we show in Figs. [9] and [10] this approximate formula does 
quite a good job of fitting the numerical results for a wide 
range of models. The only region of parameter space where 
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FIG. 9. — Numerically calculated growth rates of modes with m = 1 and 
^?sing = 1 ■ 1/f j ■ These modes have among the largest growth rates. The solid 
lines show the results for a series of models for fixed 7,,, and varying R,„/Rj. 
The dotted lines are the growth rates predicted by eq. )8U . Note the very 
good agreement except near the top of the plot, where the models are non- 
relativistic. The dashed lines are the numerical growth rates for modes with 
k = 0. 
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Fig. 10. — Numerically calculated growth rates of modes with m = 1 and 
^siiis = The solid lines show the results for a series of models with 

fixed R m/R j and varying 7,,, . The dotted lines are the growth rates predicted 
by eq. )81t . The dashed lines are the numerical growth rates for modes with 
k = 0. 

the formula fails is when the underlying equilibrium becomes 
"non-relativistic" and jj approaches unity. These models are 
near the upper end of Figs. |9]and[T0]and have extremely large 
growth rates. 

Although modes with R sm „ ~ Rj have the largest growth 
rates, these modes have relatively short wavelengths <C Rj 
along the z-axis (see eq. [80b . With such short wavelengths it 
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is not clear if these instabilities can grow to large amplitude. 
It is therefore interesting to consider modes with k — ► 0. Fig- 
ure!?] shows that Model B is stable as k — > 0, whereas Fig. |6] 
indicates that the k = mode in Model C is nearly as unstable 
as the fastest-growing mode. 

The dashed lines in Figs. l9l and [Tol show numerical results 
for the growth rates of modes with k = for various combi- 
nations of the model parameters j m and R,„. For small values 
of R m ^ 0.11?;, the results are nearly identical to those we 
described above for the fastest-growing mode (R smg = l.lRj, 
solid lines). This is to be expected based on the results shown 
in Fig. |6]for Model C, which has R m = QARj. With increasing 
R m , however, the k = modes become less unstable than the 
modes with R smg ~ Rj. By R„ 



Q3Rj, the k ■■ 



■ modes are 



fully stable, thus explaining the result shown in Model B (Fig. 
Hi, which has R m = 03Rj. 

4.5. Spatial Growth of Unstable Modes 

The discussion so far was limited to modes with real k and 
complex u>. An equally interesting problem is to consider 
modes with real ui and complex k. From equation (f37T >. we 
see that the eigenfunctions take the form 

E x = [E 1R (R)R+E 14 ,(R)$+E u (R)z] 
x exp[-i'o;/ + imcj) + iRe(k)z - lm(k)z] 
ocexp[/Re(fc)z]exp(z/Z), (82) 
where to is real and Z = -1 j\m(k) is the scale length on which 
the mode e-folds in the z-direction. Such spatially growing 
"convected" modes are particularly relevant for sources with 
long-lived steady jets. 

As discussed in iPayne & Cohnl (1 19851) and 
Appl & Camenzjnd] (119921) . there is a strong symmetry 
between modes with real k and complex u>, and those with 
complex k and real ui. In particular, the growth rates of the 
two kinds of modes are related by 

lm(k) = -lm(uj)/vg, (83) 
where v g = <9Re(w) /dRe(k) is the group velocity of the mode. 
Our unstable m = 1 modes have v g very nearly equal to c. 
Therefore, we immediately obtain from equation ( T8lT > the fol- 
lowing estimate of the spatial e-folding scale Z of the fastest- 
growing convected mode: 



Mode with R sing : 



l.IRi : — 



2.5 T/ 



(84) 



, v (l-2R m /Rj) 
Zero-frequency modes (u> = 0) should have almost the same 
Z for small values of R m /Rj, but the growth should cut off at 
a somewhat smaller value of R m compared to the modes with 
^sing = 1 • IRp 

Figure [TT] shows numerical results. Modes with ^ s i„ g = 
l.lRj have growths consistent with equation d84l i. and the 
modes with uo = have similar growths except that the insta- 
bility cuts off at somewhat smaller values of R m . The results 
are as expected and are very similar to those shown in Fig. [9] 

The above results correspond to an idealized cylindrical jet. 
In the case of real jets we must allow for a finite opening an- 
gle 9j = dRj/dz- Many force-free jet models have opening 
angles th at vary inversely as the Lorentz factor: 9j ~ few/7; 
(TMN08). Using this scaling we can estimate approximately 
the evolution of the mode amplitude a with distance: 
da 1 da 7; da 7; a 



dRj 9j dz 



few dz 



few Z 



1 



few 2.5jjRj few x 2.5 Rj 



IS] 



o 




log (R m /Rj) 

FIG. 1 1 . — Numerically calculated e-folding scale Z for modes with m = 1, 
w real and R smg = l.lRj. These modes have among the largest largest growth 
rates. The solid lines show the results for a series of models with a given 
value of 7,,, and diffe rent values of R„,/Rj. The dotted lines are the growths 
predicted by eq. j84t . Note the very good agreement except near the bottom 
of the plot, where the models are non-relativistic. The dashed lines are the 
numerical values of Z for modes with u = 0. 



Solving this differential equation, and using Rj oc z 
(eq. [6), we obtain 



a/4 , 



a{z) oc z e 



,0.5-1 



(86) 



where e is a small number < 0.1. This estimate is very crude, 
but it does suggest that, in realistic jets, the unstable kink 
mode we have studied in this paper grows only weakly with 
increasing distance. 

4.6. Towards an Improved Instability Criterion 

In §1, we introduced three different instability criteria, of 
which the IPL and TMT criteria refer specifically to rotating 
force-free jets. Since B z is practically constant in our equilib- 
ria, all of our models are close to the boundary between stabil- 
ity and instability according to the IPL criterion (eq. Sim- 
ilarly, since B 2 ^ m Vt 2 R 2 B 2 p in our equilibria, our models are 

marginally stable according to the TMT criterion (eq. [3]). 14 It 
is thus not possible to understand from either of these criteria 
why Model A is stable and Models B and C are unstable. 

It is, of course, not surprising that the above instability cri- 
teria fail. Our jet equilibria include the effects of poloidal field 
curvature, which were not considered by the previous authors. 
For easier comparison with previous work, let us rewrite our 
balance condition ( 1221 as follows: 



1 d 
WdR 



E 2 0R )R 2 



8- 



d 

dR 



-B 2 



4irR r 



(87) 

If we leave out the last term , this is equivalent to e quation (6) 
in llstomin & Parievl d 19961) and equation (16) in iLvubarskiil 
(119991) . The IPL instability criterion states that the quantity 

14 A strict application of the TMT criterion would indicate that our models 
are unstable, since > E 2 = £l 2 R 2 B 2 p . However, B 2 ^ - E 2 <§; B\, so the 
models deviate only slightly from marginal stability. 
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dB\ 7 /dR should be negative. We might wish to generalize 
this by saying that the right hand side of equation d87l i. in- 
cluding the poloidal curvature term, should be positive. Un- 
fortunately, this simple modification is not sufficient since the 
right hand side is positive for all of our models, whereas not 
all our models are unstable; not only should the right hand 
side be positive, its magnitude should be larger than some 
amount. The same seems to be true with the TMT criterion. 
This criterion indicates that all our equlibria should be unsta- 
ble, whereas only some of them are. 

Qualitatively, what distinguishes the unstable Models B and 
C from the stable Model A is that the former have made the 
transition to the second accleration regime. This is reflected 
in their 7(7?) profiles (Fig. 2) which have d 7 / dR < at larger 
radii. Thus, one might guess that instability requires the jet to 
be in the second acceleration regime and/or the jet to have a 
declining j(R). Once again, these conditions by themselves 
are not sufficient. To have an instability, -f(R) should decline 
over a sufficiently broad range of radius, e.g., R m should be 
less than ~ 0.45/?, in our models. 

An alternate approach which we have found useful is to 
focus on the left-hand side of equation d8Tb . From equa- 
tions ( TTol l. ( TP7I ). ( fT8l , we see that for our equilibria we have 
Bq^-Eqx = (R/R m ) 4 BQ„. Furthermore, we have seen that 
modes with R s \ ng ~ Rj (the fastest-growing modes) are unsta- 
ble so long as R,„ < 0A5Rj, while modes with k — > (long- 
wavelength modes) are unstable for R m < 03Rj. From this, 
we obtain the following approximate instability criteria: 

Modes with R sing ~ Rj : (B^-E^) V * > 5|B 0z |, (88) 
Modes with k -> : (B 2 0ct> -E? )R ) 1/2 > 12|B 0z |, (89) 

where we have set R = Rj to obtain the numerical coefficients 
on the right. These conditions are easier to interpret if we 
boost to the comoving frame of the fluid (V. Pariev, private 
communication), where the electric field vanishes. In this 
frame, the criterion for instability becomes 

Modes with R sing ~ Rj : |B O 0,comov I > 5 \B Qz . mmm | , (90) 
Modes with k — > : |Bo0,comov| > 12|B (te ,comov|- (91) 

That is, in the comoving frame, the toroidal field must domi- 
nate the poloidal field by more than a certain critical factor. 15 
Written in this form, the condition resembles the KS criterion 
(eq.[B. 

5. SUMMARY AND DISCUSSION 

The relativistic jet model we have considered in this paper 
is particularly simple: it is cylindrical, it assumes force-free 
conditions, and it assumes rigid rotation. Within the limita- 
tions of these reasonable approximations, we have attempted 
to be as close to numerically simulated jets as possible. We 
include the effect of poloidal field curvature, which is known 
to play an important role in numerical force-free jets (§2.1), 
and we choose functional forms for the various field compo- 
nents in the equilibrium (§2.2) to m atch as clo sely as possible 
our previous force-free simulations (TMNOjl). 

Our equilibrium model is described by two parameters: the 
maximum Lorentz factor j m , and the radius at which this max- 
imum is achieved R m . The ratio of the latter to the jet radius 

' 5 Sin ce o ur equilibria assume a constant f2 for all field lines, the criteria 
{90} and )9U are technically valid only for such models. However, since the 
criteria have been written without any explicit reference to Q, they may be 
valid more generally even when f2 varies with R. 



Rj determines the basic physics of the equilibrium. Models 
in which R m /Rj > 1 have 7(7?) increasing monotonically with 
radius R out to some maximum Lorentz factor jj < J m at the 
outer edge of the jet. Model A (Fig. [2]) is an example. In 
these models the e ntire jet i s in the first acceleration regime 
(see §2.1,2.2.1 and TMN08 for details). We find that all these 
models are perfectly stable. 

Models with R m /Rj < 1 are more interesting. Here, j{R) 
increases upto a maximum value j m at R = R,„ and then de- 
creases down to a Lorentz factor 7/ < j m at R = Rj. Mod- 
els B and C (Fig. |2|l are examples of this kind of model. In 
these models, the jet fluid at R < R m is in the first accelera- 
tion regime, while the fluid at R,„ < R < Rj is in the second 
acceleration regime. We find that the subset of these models 
with R m /Rj < 0.45 are linearly unstable. For R,„/Rj just be- 
low 0.45, all the unstable modes have short wavelengths in 
the z-direction: A = 2?r/k z ~ 2nRj/jj. With decreasing R,„, a 
wider range of k z becomes unstable, and for R m /Rj < 0.3, we 
find that waves with k z = 0, i.e., with arbitrarily long wave- 
lengths, are unstable. The latter modes are perhaps of most 
interest since they are likely to grow to the largest amplitudes. 
The numerical results are summarized in Figs. [314TTI 

The unstable modes we find are all kink modes with az- 
imuthal wavenumber m = ± 1 . These are non-axisymmetric 
modes in which the jet is distorted helically. A key feature 
is that, at each z, the center of mass of the jet is shifted 
away from R = 0. It is well-known that MHD configurations 
with toroidal fie lds are especially susceptible to the kink mode 
(Bateman 1978), and our models follow this trend. However, 
because our equilibria both rotate and move relativistically 
along z, the criterion for instability is different from the usual 
KS criterion (eq. [T). 

We find that the typical growth rate of the unstable kink 
mode in our jet models is given by equation (T8TT >: the e- 
folding time is of order 7/ times the light-crossing time Rj/c 
across the jet. For convected modes with a real frequency, this 
translates to an e-folding length scale of order 7/ times the jet 
radius Rj. Since jets typically have opening angles ~ few/7/, 
the net result is that the unstable modes grow only slowly with 
distance from the base of the jet (eq. [86b . Of course, relativis- 
tic jets in astrophysical sources propagate over many decades, 
so in principle even this slow growth might lead to a large 
amplitude of the perturbation. Nevertheless, the fact that the 
growth is very slow reduces the seriousness of the kink insta- 
bility. 

Our jet equilibria turn out to be close to the boundary be- 
tween stability and instability according to either the IPL or 
TMT criterion (eqs. [2] [3j, so these criteria are not useful for 
interpreting the results. In addition, since our models include 
the effects of poloidal field curvature, they lie outside the 
range of validity of the IPL and TMT criteria. The most use- 
ful instability criterion we have come up with is that, in the 
comoving frame of the jet fluid, the tangential field should be 
an order of magnitude or more larger than the poloidal field 
(eqs. l90ll9Tl ). Expressed thus, the criterion is similar to the KS 
criterion (Q]), except that it should be applied in the comoving 
frame and z should be taken to be ~ Rj. 

All the work described here assumes a rigid wall enclosing 
the jet at the boundary R = Rj. We have done some calcula- 
tions with a constant pressure boundary and we find unstable 
modes with much larger growth rates compared to the rigid 
wall case. However, since we are dealing with a force-free 
jet, it is not clear that a constant pressure boundary is particu- 
larly meaningful. For instance, if the pressure is from a non- 
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relativistic gaseous envelope or cocoon, the gas would have 
substantial inertia and (we suspect) would probably behave to 
first approximation like a rigid wall. 

Various authors have discussed mechanisms by which 
instabilities might be suppressed in astrophysical jets. 
lHardee et al.1 (|2007, and references therein) have shown 
that an external wind or cocoon can stabilize the Kelvin- 
Helmholtz mode in MHD jet s, though it is not clear if this 
is rele vant for force-free jets. Mo ll. Spruit. & Obergaulinger 
(2008) show that lateral expansion causes instabilities to grow 
more slowly. In a sense, we have already included this effect 
when we derived the growth rate estimate given in equation 
(f86b . In ad dition, we n ote that some of the growth suppres- 
sion seen by Moll et al. is probably because expansion causes 
different parts of the jet to lose causal contact with one other. 
This is not an issue for force-free models, where signals prop- 
agate at the speed of light. 

It would be interesting to simulate numerically the unstable 
modes described in this paper. Apart from verifying the linear 
theory, such calculations will reveal the non-linear develop- 
ment of the mode. Does the kink mode saturate at a finite 



amplitude and lead to a more-or-less coherent helical pattern 
or does it destroy the initial equilibrium? This important ques- 
tion can be answered only with 3D simulations. Since the kink 
mode involves lateral motion of the jet across the axis R = 0, 
the numerical technique used must be flexible enough to allow 
such motions (e.g., as described by iMcKinn ev & Blandford 
l2008h . 

We conclude by reminding the reader that the work de- 
scribed here refers to a particularly simple model of relativis- 
tic jets which is based on the force-free approximation. In real 
jets, once the flow crosses the fast magnetosonic point, the in- 
ertia of the gas starts to play a role and the force-free appr oxi- 
mation is no longer valid (e.g.. lTchekhovskoy et alJl2009l) . In 
this regime, we must consider the full MHD equations. 
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APPENDIX 

A. CONSTANCY OF POLOIDAL MAGNETIC FIELD ACROSS FORCE-FREE JETS 

Figure [2j] shows that in numerical simulations of force-free jets B p hardly changes with R. In this Appendix we show that this 
is a common feature of all jet solutions that smoothly connect to a spinning compact object at the base. 

Consider the force balance equation (0. Sufficiently near the compact object, where the jet is in the first acceleration regime 
(see 32.11 ). we can drop the terms proportional to R~ l and (B\— E 2 ) since in the first acceleration regime j 2 <C y 2 leading to 

E 2 /B 2 p <C R c /R and B^-E 2 < B 2 dTMN08l) . Then the force balance equation © simplifies to 



d(B 2 p ) 
dR 



0. 



(Al) 



Therefore, sufficiently near the compact object the poloidal field is nearly constant, 

B p (R) w const. (A2) 
Each field line is labeled by the amount of poloidal magnetic flux <f> that it encloses. Due to ( IA2I ) this flux can be written simply 



as 



3> 



ttB p R 2 . 



(A3) 



These relations are valid throughout the first acceleration regime. We now show that they actually hold asymptotically in all parts 

of the j et. 

Recall that in force-free magnetosphere s the enclosed poloi dal c urrent / is preserved along each field line dMestell fl 96 lh 
Okam otolll978HThorne et aiTl 986; Beskrnl ll997l:lNaravan et al1l2007l: |Tchekhovskoy et al. 2008), 



/ = /($)= -RB 4 



■B n R 



(A4) 



where the approximate equality is due to B^ « —E ■ 
that 57 is conserved along field lines, we obtain 



-flRB p /c for R^> Ra (see eq.fTTTi. Comparing ( IA3b and ( lA4b and recalling 



2tt 



(A5) 



Since both sides of this relation depend onl y on <£>, this relation is valid eve rywh ere in the s olutio n, e ven th ough we derived it only 
in the first acceleration regime. Using (IA4b to substitute for / yields back dA3t . Thus eqs. (IA2b and ( IA3b are valid everywhere in 
the jet. 

B. COEFFICIENTS IN EQUATIONS (52) AND (53) 
In this Appendix we give explicit expressions for the coefficients D\(R)-D(,(R) defined in §3.1. The functions are 

cm 



D\(R) = 
D 2 (R) = 



ck 
— ! 



luR 



cm cm f, cm 
ujR z ujR ujR 



, ick icm iu> 
D 3 (R)=[ + ^ 

uj ujR 1 c 



(Bl) 
(B2) 

(B3) 
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ck cm D 3 
D 4 (R)=-f+—g-h-C 3 -^, 
lu luR D\ 



C 5 



D 5 (R) = C 4 + C 3 -^+[C b + ^-)c 1 +C 5 C' 1 -C 7 -^- 



R 



D 2 



+C 3 



D 2 D[ D> 2 \ ic 
D\ Dj w 8 



D 2 D 3 D 3 (D 3 D\ D\ 

D6{ r )=Cs+ c 3 ^-c 71 ± + c 3 ^i-- 



where primes denote derivatives with respect to R, and the functions C\ (R)-Cn(R) are given by 

ujRg - cmh 



Ci(R) = 
C 2 (R) = 



uiRf-ckRh' 
ck cm 

Zr^'^r 2 ^ 



h h ck , , 

+ —f-h, 

R R r lu 



C 3 (R)=-(f+gC l ), 

LU 



C 4 (R) = 
C 5 (R) = 

C 6 (R) = 
Ci{R) = 
Cs(R) = 



icm 2 iu> ickm 

- f+ - f+ UR S - ikh > 



' luR 2 ' 



IC IC 
LuR 8 UjRr 



IC 



lu 



ickm 



ick 2 



2ic iuj im 
ujR z c R 



luR 

2ic ic . 2ic . 

■J+-f' + C 5 C l + —gC' 1 



ujR 
2ck 



LU 

cm 



cm 2h 
luR-'~Zr 2 ~ 8 ~~R + ujR c R 



h ck , cm , , 
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Rr LU LUR 



(B4) 

(B5) 
(B6) 

(B7) 
(B8) 
(B9) 
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(Bll) 
(B12) 
(B13) 
(B14) 
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